A Quick Start Tutorial¶
This tutorial offers a quick-start guide for using BaySpec to fit a spectral model to gamma-ray data. It can be broadly divided into the following three sections:
Data: The spectra from Fermi/GBM’s NaI detector and BGO detector.
Model: A simple cutoff power-law function.
Fitting: Bayesian inference implemented using multinest.
import numpy as np
from bayspec.model.local import *
from bayspec import DataUnit, Data, Infer, Plot
Load spectra and response data.
nai = DataUnit(
src='./ME/me.src',
bkg='./ME/me.bkg',
rsp='./ME/me.rsp',
notc=[8, 900],
stat='pgstat',
grpg={'min_sigma': 3, 'max_bin': 10})
bgo = DataUnit(
src='./HE/he.src',
bkg='./HE/he.bkg',
rsp='./HE/he.rsp',
notc=[300, 38000],
stat='pgstat',
grpg={'min_sigma': 3, 'max_bin': 10})
data = Data([('nai', nai),
('bgo', bgo)])
2.Define spectral model.
model = cpl()
3.Run bayesian infer.
infer = Infer([(data, model)])
print(infer)
╒════════╤══════════════╤═════════════╤═════════════╤═════════╕
│ cfg# │ Expression │ Component │ Parameter │ Value │
╞════════╪══════════════╪═════════════╪═════════════╪═════════╡
│ 1 │ cpl │ cpl │ redshift │ 0 │
╘════════╧══════════════╧═════════════╧═════════════╧═════════╛
╒════════╤══════════════╤═════════════╤═════════════╤═════════╤═════════════╕
│ par# │ Expression │ Component │ Parameter │ Value │ Prior │
╞════════╪══════════════╪═════════════╪═════════════╪═════════╪═════════════╡
│ 1* │ cpl │ cpl │ $\alpha$ │ -1 │ unif(-8, 4) │
├────────┼──────────────┼─────────────┼─────────────┼─────────┼─────────────┤
│ 2* │ cpl │ cpl │ log$E_{c}$ │ 2 │ unif(0, 4) │
├────────┼──────────────┼─────────────┼─────────────┼─────────┼─────────────┤
│ 3* │ cpl │ cpl │ log$A$ │ -1 │ unif(-6, 5) │
╘════════╧══════════════╧═════════════╧═════════════╧═════════╧═════════════╛
post = infer.multinest(nlive=300, resume=True, savepath='./multinest')
*****************************************************
MultiNest v3.10
Copyright Farhan Feroz & Mike Hobson
Release Jul 2015
no. of live points = 300
dimensionality = 3
resuming from previous job
*****************************************************
Starting MultiNest
Acceptance Rate: 0.689121
Replacements: 6271
Total Samples: 9100
Nested Sampling ln(Z): -229.513355
Importance Nested Sampling ln(Z): -229.263794 +/- 0.024273
ln(ev)= -229.21088552733494 +/- 0.24113203804165545
Total Likelihood Evaluations: 9100
Sampling finished. Exiting MultiNest
analysing data from ./multinest/1-.txt
print(post)
╒════════╤══════════════╤═════════════╤═════════════╤════════╤══════════╤════════╤══════════════════╕
│ par# │ Expression │ Component │ Parameter │ Mean │ Median │ Best │ 1sigma CI │
╞════════╪══════════════╪═════════════╪═════════════╪════════╪══════════╪════════╪══════════════════╡
│ 1 │ cpl │ cpl │ $\alpha$ │ -1.563 │ -1.563 │ -1.561 │ [-1.572, -1.552] │
├────────┼──────────────┼─────────────┼─────────────┼────────┼──────────┼────────┼──────────────────┤
│ 2 │ cpl │ cpl │ log$E_{c}$ │ 2.69 │ 2.69 │ 2.688 │ [2.673, 2.707] │
├────────┼──────────────┼─────────────┼─────────────┼────────┼──────────┼────────┼──────────────────┤
│ 3 │ cpl │ cpl │ log$A$ │ -0.771 │ -0.771 │ -0.77 │ [-0.778, -0.765] │
╘════════╧══════════════╧═════════════╧═════════════╧════════╧══════════╧════════╧══════════════════╛
╒════════╤═════════╤═════════════╤════════════╤════════╕
│ Data │ Model │ Statistic │ Value │ Bins │
╞════════╪═════════╪═════════════╪════════════╪════════╡
│ nai │ cpl │ pgstat │ 388.00 │ 106 │
├────────┼─────────┼─────────────┼────────────┼────────┤
│ bgo │ cpl │ pgstat │ 32.62 │ 26 │
├────────┼─────────┼─────────────┼────────────┼────────┤
│ Total │ Total │ stat/dof │ 420.62/129 │ 132 │
╘════════╧═════════╧═════════════╧════════════╧════════╛
╒════════╤════════╤════════╤═════════╕
│ AIC │ AICc │ BIC │ lnZ │
╞════════╪════════╪════════╪═════════╡
│ 426.62 │ 426.81 │ 435.27 │ -229.26 │
╘════════╧════════╧════════╧═════════╛
fig = Plot.infer_ctsspec(post, style='CE', show=True)
fig = Plot.post_corner(post, show=True)
earr = np.logspace(np.log10(0.5), 3, 100)
modelplot = Plot.model(ploter='plotly', style='vFv', CI=True)
fig = modelplot.add_model(model, E=earr, show=True)